//
// Ephi - simulation of magnetic fields and particles
// Copyright (C) 2007 Indrek Mandre <indrek(at)mare.ee>
// For more information please see http://www.mare.ee/indrek/
//
// This program is free software; you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation; either version 2 of the License, or
// (at your option) any later version.
//
// This program is distributed in the hope that it will be useful,
// but WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
// GNU General Public License for more details.
//
// You should have received a copy of the GNU General Public License along
// with this program; if not, write to the Free Software Foundation, Inc.,
// 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
//

#include <math.h>
#include <assert.h>
//#include <sys/time.h>

#include "ephi.hpp"

LinearSystem::LinearSystem (size_t n, size_t k) : n(n), k(k), npk(n+k)
{
  rows.resize (n);
  data = new prec_t [n * npk];
  for ( size_t i = 0; i < n; i++ )
      rows[i] = data + i * npk;
  results = new prec_t [n * k];
}

LinearSystem::~LinearSystem ()
{
  delete [] data;
  delete [] results;
}

namespace
{
  struct gjtask : Task
  {
    LinearSystem *ls;
    size_t npk;
    prec_t **begin;
    prec_t **end;
    prec_t **except;
    prec_t *fset;
    size_t i;
    // this is not as good as it seems, the math here is so trivial it takes almost no time at all
    // as the results most of the time is done waiting behind the memory - and that is slow.
    void execute ()
    {
      for ( prec_t **rit = begin; rit != end; rit++ )
        {
          if ( rit == except )
              continue;

          prec_t *f = fset;
          prec_t *it = *rit + i;
          prec_t *end = *rit + npk;
          prec_t mod = *it;
          *it = 0;
          it++;

          while ( it != end )
            {
              *it -= mod * *f;
              it++;
              f++;
            }
        }
    }
  };

  struct gbstask : Task
  {
    size_t npk;
    prec_t **begin;
    prec_t **end;
    prec_t **except;
    prec_t *fset;
    size_t i;
    // this is not as good as it seems, the math here is so trivial it takes almost no time at all
    // as the results most of the time is done waiting behind the memory - and that is slow.
    void execute ()
    {
      for ( prec_t **rit = begin; rit != end; rit++ )
        {
          prec_t *f = fset;
          prec_t *it = *rit + i;
          prec_t *end = *rit + npk;
          prec_t mod = *it / *fset;
          *it = 0;
          it++;
          f++;

          while ( it != end )
            {
              *it -= mod * *f;
              it++;
              f++;
            }
        }
    }
  };
};

bool LinearSystem::solveGJ ()
{
  size_t tc = THREADS_CPUCOUNT;
  TaskManager tmgr(tc);
  std::vector<gjtask> tasks;
  tasks.resize (tc);
  size_t pc = n / tc;
  if ( !pc ) pc = 1;
  for ( size_t i = 0; i < tasks.size(); i++ )
    {
      tasks[i].ls = this;
      tasks[i].npk = npk;
      tasks[i].begin = getRows() + i * pc;
      tasks[i].end = tasks[i].begin + pc;
      if ( tasks[i].begin > getRows() + n )
          tasks[i].begin = getRows() + n;
      if ( i == tasks.size() - 1 || tasks[i].end > getRows() + n )
          tasks[i].end = getRows() + n;
    }

  for ( size_t i = 0; i < n; i++ )
    {
      //printf ("LinearSystem:: solveGJ %u/%u\n", (unsigned int)i, (unsigned int)n);
      // find the greatest pivot and switch rows
      size_t m = i;
      for ( size_t j = i + 1; j < n; j++ )
        {
          if ( fabs(rows[j][i]) > fabs(rows[m][i]) )
              m = j;
        }
      if ( m != i )
        {
          prec_t *r = rows[m];
          rows[m] = rows[i];
          rows[i] = r;
        }

      //printf ("pivot:\n");
      //print();

      // now divide the current row by first nonzero element
      prec_t *it = rows[i] + i;
      prec_t *end = rows[i] + npk;
      prec_t div = *it;
      if ( fabs(div) < 1e-23 )
          return false;
      *it = 1.0;
      it++;
      prec_t *fset = it;
      //struct timeval tv1, tv2, tv3;
      //gettimeofday (&tv1, 0);
      while ( it != end )
        {
          *it /= div;
          it++;
        }
      //gettimeofday (&tv2, 0);

#if 1
      for ( size_t j = 0; j < tasks.size(); j++ )
        {
          tasks[j].except = &rows[i];
          tasks[j].i = i;
          tasks[j].fset = fset;
          tmgr.addTask (&tasks[j]);
        }

      while (tmgr.run());
#else
      for ( size_t j = 0; j < n; j++ )
        {
          if ( j == i )
              continue;

          prec_t *f = fset;
          prec_t *it = rows[j] + i;
          prec_t *end = rows[j] + npk;
          prec_t mod = *it;
          *it = 0;
          it++;

          while ( it != end )
            {
              *it -= mod * *f;
              it++;
              f++;
            }
        }
#endif
      //gettimeofday (&tv3, 0);
      //printf ("%zu, %zu\n", (size_t)(tv2.tv_usec - tv1.tv_usec) + (size_t)(tv2.tv_sec - tv1.tv_sec)*1000000, size_t(tv3.tv_usec - tv2.tv_usec) + (size_t)(tv3.tv_sec - tv2.tv_sec)*1000000);
      //printf ("result:\n");
      //print();
    }

  for ( size_t j = 0; j < k; j++ )
    {
      prec_t *rb = results + j * n;
      for ( size_t i = 0; i < n; i++ )
        {
          *rb++ = rows[i][n + j];
        }
    }

  return true;
}

bool LinearSystem::solveGBS ()
{
  size_t tc = THREADS_CPUCOUNT;
  TaskManager tmgr(tc);
  std::vector<gbstask> tasks;
  tasks.resize (tc);
  for ( size_t i = 0; i < tasks.size(); i++ )
      tasks[i].npk = npk;
  for ( size_t i = 0; i < n - 1; i++ )
    {
      Ephi::debug (Ephi::DEBUG_INFO, "LinearSystem:: solveGBS %u/%u\n", (unsigned int)i + 1, (unsigned int)n);
      // find the greatest pivot and switch rows
      size_t m = i;
      for ( size_t j = i + 1; j < n; j++ )
        {
          if ( fabs(rows[j][i]) > fabs(rows[m][i]) )
              m = j;
        }
      if ( m != i )
        {
          prec_t *r = rows[m];
          rows[m] = rows[i];
          rows[i] = r;
        }

      //print();

      prec_t *fset = rows[i] + i;
#if 1
      // divide the calculation space across threads
      prec_t **rbegin = &rows[0] + i + 1;
      prec_t **rend = &rows[0] + n;
      size_t pc = (rend - rbegin) / tc;
      if ( !pc ) pc = 1;
      for ( size_t j = 0; j < tasks.size(); j++ )
        {
          tasks[j].fset = fset;
          tasks[j].i = i;
          tasks[j].begin = rbegin + j * pc;
          tasks[j].end = tasks[j].begin + pc;
          if ( tasks[j].begin > rend )
              tasks[j].begin = rend;
          if ( j == tasks.size() - 1 || tasks[j].end > rend )
              tasks[j].end = rend;
          tmgr.addTask (&tasks[j]);
        }
      while (tmgr.run());
#else
      for ( size_t j = i + 1; j < n; j++ )
        {
          prec_t *f = fset;
          prec_t *it = rows[j] + i;
          prec_t *end = rows[j] + npk;
          prec_t mod = *it / *fset;
          *it = 0;
          it++;
          f++;

          while ( it != end )
            {
              *it -= mod * *f;
              it++;
              f++;
            }
        }
#endif
      //printf ("step:\n");
      //print();
    }

  // and now solve the results
  for ( size_t j = 0; j < k; j++ )
    {
      prec_t *rb = results + j * n + n - 1;
      prec_t **rit = &rows[n-1];
      prec_t **rend = &rows[0] - 1;
      size_t l = n - 1;
      while ( rit != rend )
        {
          prec_t *row = *rit;
          prec_t *it = row + l;
          prec_t *end = row + n;
          prec_t cum = *(end + j);
          prec_t div = *it++;
          prec_t *tail = rb + 1;
          while ( it != end )
            {
              cum -= *it * *tail;
              it++;
              tail++;
            }
          rit--;
          *rb-- = cum / div;
          l--;
        }
    }

  return true;
}

void LinearSystem::print ()
{
  const char *fmt1;
  if ( n <= 10 )
      fmt1 = "%u:";
  else if ( n > 10 )
      fmt1 = "%2u:";
  else if ( n > 100 )
      fmt1 = "%3u:";
  else if ( n > 1000 )
      fmt1 = "%4u:";
  else if ( n > 10000 )
      fmt1 = "%5u:";
  else
      fmt1 = "%8u:";
      
  for ( size_t i = 0; i < n; i++ )
    {
      printf (fmt1, (unsigned int)i);
      for ( size_t j = 0; j < npk; j++ )
          printf (" %.2e", prec2double(rows[i][j]));
      printf ("\n");
    }
}

void LinearSystem::clear ()
{
  memset (data, 0, sizeof (prec_t) * n * npk);
}

